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ABSTRACT 

We calculate the global quasi-steady state of a thin disk perturbed by a low-mass protoplanet 
orbiting at a fixed radius using extremely high-resolution numerical integrations of Euler's equations 
in two dimensions. The calculations are carried out using a moving computational domain, which 
greatly reduces advection errors and allows for much longer time-steps than a fixed grid. We calculate 
the angular momentum flux and the torque density as a function of radius and compare them with 
analytical predictions. We discuss the quasi-steady state after 100 orbits and the prospects for gap 
formation by low mass planets. 

Subject headings: hydrodynamics - methods: numerical - planet-disk ineractions - planets and sat- 
tellites: formation - protoplanetary disks 
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1. INTRODUCTION 

It has been well over a decade since the discovery of 
"Hot Jupiters" , i.e. extrasolar planets orbiting their host 
stars at radii at which they are not expected to be able 
to form ([Mayor fc Oueloz 1119951 ). Yet the study of pro- 
toplanetary migration still poses fundamental theoreti- 
cal challenges. Migration is expected to be driven by 
gravitational to rques exerted on the pla net by the proto- 
planetary disk (jGoldreich fc Tremaine IIT98 01. A major 
unsolved puzzle is that expected migration timescales are 
short compared to the lifetime of the disk (Type I Migra- 
tion) (Goldreich fc Tremaine ' 1979, 1980; Ward 'l997j; 
Armitage fc Rice .2005,; Papalpizou fc Terquem 20061 ; 
Yu et al. Il2010[ ). This suggests that a mechanism to halt 



planetary migration (or at least slow it down) is required, 
otherwise planets are expected to migrate into their pro- 
tostar or be kicked out of the system. 

If a planet opens a gap in the disk, however, migration 
slows down to viscous timescales (T ype II Migration) 
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Crida ct al. 2007). Depending on disk viscosity, this 
could significantly increase the time it takes for the 
planet to migrate. It is therefore of interest to precisely 
determine what conditions are required for a planet to 
open a gap. 

The question of criteria f or gap formation h a s bee n 
explored previously, e.g . 'Ward fc Hourigan] ([1989[ ): 
iLin fc Papaloizou I (|l993l) : [Rafikov (2002b). A gap is 
formed when gas is driven away from the planet, which 
requires an exchange of angular momentum between the 
planet and the gas in its neighborhood. Merely exuding 
torque in the form of a spiral density wave (Figure [1]) is 
not sufficient - this torque must actually be deposited 
in neighboring fluid elements, thus transporting angular 
momentum and driving disk evolution. In the limit of 
a very low-mass perturber, the density wave driven by 
the planet will simply provide a uniform flux of angular 
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Fig. 1. — The perturbation to surface density caused by a low- 
mass planet, Mp = .0209Afrfi. The planetary wake traces out a 
spiral shape. This calculation used 4096 radial zones. 

momentum, without opening a gap. 

On the other hand, if the perturbing mass is large 
enough to create a strongly nonlinear wave, the wave will 
quickly shock, dissipate, and transfer angular momentum 
between the planet and the gas, driving material away 
from the planet and forming a gap. Between the two ex- 
tremes there must exist a threshold mass for which a gap 
forms. The threshold mass, Mth, can be estimated by 
the condition for strong nonlinearity in the spiral wave 
([Lin fc Papaloizou |[T993[) : 



(1) 



where c is the sound speed, J7 is the orbital frequency, 
r is the orbital radius, and M. = rpftp/cp is the Mach 
number, which for a thin disk is assumed to be equal 
to r/h, where h is the disk scale height (the subscript 
"p" means the quantities are evaluated in the vicinity 
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of the planet). This threshold mass is also the mass 
for which the planet's Hill radius is of order the disk 
scale height. If we assume the central star is a solar 
mass and the disk is orbiting at Mach 20, the thresh- 
old mass is 41Af^. Atlernatively, assuming a minimum 
mass solar nebula, the threshold mass is approximately 
21Me(r/lAf7)3/4 (|HavashillT98l . In our work we use 
the former value (the two answers agree at = 2.4At7). 
This estimate gives an upper bound; strong nonlinearity 
leads to gap opening. However, the actual gap-opening 
threshold might actually be smaller than this estimate. 
The basic idea behind ([T]) is correct. Gap opening re- 
quires not only the generation of a spiral wave, but a 
damping mechanism which transfers angular momentum 
from the wave to local fluid elements. However, a point 
raised by Goodman & Rafikov (2001) is that there is 
damping even in the weakly nonlinear regime. In other 
words, ([1]) is a sufficient condition, but it might not be a 
necessary one. A weakly nonlinear perturbation will also 
shock eventually, some distance from the planet. Assum- 
ing negligible disk viscosity and that the wave shocks 
before reaching the edge of the disk, any massive per- 
turber at a fixed radius is capable eventually forming a 
gap, though for small enough perturbers the timescale is 
prohibitively long and Type I migration would occur too 
rapidly for gap opening. Hence, for low mass planets, the 
formation of a gap depends on a competition between the 
timescale for gap opening and other relevant timescales, 
such as the viscous timescale or migration timescale. 

In order to examine the details of ga p opening by a 
weak shock. [Goodman fc Rafikovl (|2001l ) determined the 
wave form produced by a planet in the linear regime, then 
described its evolution into a shock assuming weak non- 
linearity. This semi-analytic result was calculated in the 
shearing box approximation, and then later too k into ac- 
count the cylindrical geometry of the system ([Rafikov I 
I2002a|) . For the latter case, some stronger predictions 
were made; in the global analysis, it was possible to take 
into account large-scale variations in density and sound 
speed in the disk, and to determine the conditions un- 
der which the density wave would leave the disk without 
shocking. 

Recently, numerical investigations have confirmed the 
semi- analytic pred i ctions in the shearing box approxima- 
tion. iMuto et al. I ([2010( ) reported a dip in d ensity ("par- 
tial gap " ) around planets of mass ~ .2MTh- [Dong et al. I 
([2011aD performed a higher resolution calculation, fo- 
cusing on smaller planetary masses and recovering the 
semi-analytical predictions to high accuracy. Until now, 
however, no global calculations have been performed to 
study the disk behavior, in particular directly calculat- 
ing gap opening by this mechanism. The global case 
is much more challenging numerically, for several rea- 
sons. First, there is simply a much larger computa- 
tional domain, which is not likely to be circumvented 
by using higher resolution near the planetary orbit, if 
the wave shocks far from the planet. Secondly, the rele- 
vant dynamical timescales are sound-crossing timescales 
('^ r/cs), but because protoplanetary disks orbit super- 
sonically, the time-steps are generally Courant-limited 
by the orbital timescales at the innermost resolved orbit. 
This can be orders of magnitude shorter, which means it 
may require a prohibitively large number of timesteps to 



reach a state resembling quasi-equilibrium. Thirdly, and 
perhaps most importantly, because the disk is supersonic 
the bulk of the motion is pure advection across the grid, 
and underresolved advection errors can completely wash 
out subtle features of the motion (like a weakly nonlinear 
shock forming) taking place in the local Keplerian frame. 

For these reasons, we perform our own calculations 
of proto-planetary disks using a new numerical method, 
whereby instead of using a fixed numerical grid, we al- 
low the computational cells to move and shear past 
one another with the bulk Keplerian flow. The nu- 
merical method we use is a variant of the TESS code 
([Dufii'ell fc MacFadven I [2011[ ). with several important 
modifications specifically designed for disk problems. 
TESS uses moving finite volumes to solve the equations 
of gas dynamics in conservation form. The motion of 
the cells is accomplished by performing a Voronoi tes- 
sellation of the computational domain each time-step. 
However, the numerical scheme is completely specified 
for any kind of tessellation, so in principle the domain 
can be decomposed into whatever cell shapes are most 
advantageous. In the present work, we choose to de- 
compose the domain into wedge-like annular segments, 
as typically implemented for cylindrical (r, cj)) grids (see 
Figured]). The cells remain at fixed radii and rotate with 
the local angular velocity of the fluid. As a result, the 
global calculations are effectively computed on a locally 
co-moving numerical mesh. 

After reviewing relevant theoretical predictions from 
the literature in ^\ and describing pertinent details of 
our numerical techniques in fJJl we present the results 
of our calculation in 21 including a demonstration of 
convergence and a calculation of the global distribution 
of torque density and angular momentum flux, before 
summarizing in fjsj 

2. SEMI-ANALYTIC PREDICTIONS 

Here we briefly summarize the theory regarding 
the generation of the spiral density wave and its 
subsequent nonlinear evol u tion. For det ails, see 
IGoodman fc Rafikov I ([200l . IRafikov I (l2002aD . We al- 
ways work in the thin-disk approximation, where we ig- 
nore all vertically propagating modes, and the equations 
of motion reduce to Euler's equations in two dimensions: 

9tS + a,(v,I]) =0, (2) 

dt i^v, ) + d,{J:v,v, + PS,, ) = , (3) 

dt{E) + d,ME + P)) ^ F,v,. (4) 

In the above, S is the surface density, v is the fluid ve- 
locity, P is pressure, F is the total gravitational force in- 
cluding the central mass of the protostar plus the orbiting 
planetary potential and E is the fluid energy density, 

+ (5) 

We use e to denote the internal energy density of the 
fluid, which is related to the pressure via the equation of 
state, 

P = (7 - l)e, (6) 

and 7 is the adiabatic index of the fluid. In this work, we 
set 7 = 1.001, effectively producing an isothermal equa- 
tion of state. 
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2.1. Constant Density and Pressure 

The formulae for the semi-analytic predictions simplify 
if we assume a uniform background surface density and 
sound speed, and Keplerian velocity, 

So (''), Co (r) = constant, (7) 

no{r) = 17p(rp/r)3/2. (8) 

In this case, the linearized equations (jRafikov I I2002af) 
predict that the spiral density wave generated by the 
planet will trace out a path described by: 



sign(r - rp) 3 - 2 



M. (9) 



The perturbation can most conveniently be described us- 
ing the following variables (jRafikov Il2002al ): 



77 = (^0 + sign(r - rp) (^3 - 2^- ^ ) X ) , 



X = 



7 + l(5S 



1/2 



2 So I A^|(rp/r)3/2-l| 



t 



25/4^ 



r/rp 



o3/2 



(10) 
(11) 

(12) 



r] acts like an azimuthal coordinate, x acts as a proxy 
for the density perturbation ST,, and t is essentially a ra- 
dial coordinate, describing distance from the planet. The 
integral in the formula for t can be evaluated in terms of 
hypergeometric functions. When the density perturba- 
tion is expressed in terms of these variables, the weakly 
nonlinear description of the system is governed by the in- 
viscid Burgers' equation. This is true both in the global 
case, and in the shearing box approximation. Therefore, 
when described in these coordinates, the shock profiles 
computed globally should look identical to the shear- 
ing box profiles calculated in Dong. R afikov &: Stone 1 

poTTbl) . 

By assuming weak nonlinearity, Euler's equation for 
the density can be expressed in terms of rj, x, and t, 
and the resulting form ulae reduce to Burgers' equation. 
IGoodman fc Rafikov I ([2001) found an approximate for- 
mula for the location of the shock in terms of t, 



Ml 



(13) 



tsh = 1.89-1-0.79 

where Mi = {2/3)MTh- 

2.2. Angular Momentum Flux 

Planets generate an effective viscosity in the disk, 
which can be characterized by the angular momentum 
flux 

/.27r 

(14) 



Pji"^) = / {TSv^r)vrrd4> 
Jo 

This can in turn be expressed in terms of a dimensionless 
function, $(i): 




Fig. 2. — The optional non-voronoi tessellation employed by Tess 
when solving problems involving gaseous disks. Finite volumes are 
wedge-like annular segments which rotate independently at fixed 
radii. 



(16) 



Because xiv^ satisfies Burger's equation, is con- 
served until a shock forms. After the shock has com- 
pletely formed and its sha.pe is described by an "N- 
wave" ([Landau fc Lifshitz II1959I ). the integral falls off 
like t~'^^^. Therefore, a reasonable approximate formula 
for the angular momentum flux is: 



m 



t < tsh 

t > tsh 



(17) 



In reality there is some smooth transition region, after 
the initial shocking time but before the shock has com- 
pletely been converted into an N-wave. The direct solu - 
tion of Burger's equation bv lGoodman k, Rafikov I ([2001[ ) 
showed that the transition region adds a small correction 
to the appoximate scaling relation. 

For radii at which the angular momentum flux (|17|) 
is not uniform, angular momentum is being transferred 
from the density wave to the disk, which creates an effec- 
tive viscosity. Wherever the slope of (|17p is large, there 
will be disk evolution toward opening a gap. 

3. NUMERICAL METHOD 
3.1. TESS code 

TESS is a finite volume hydrodynamics code con- 
structed with a moving numerical mesh which can be- 
come effectively Lagrangian when the motion of its nu- 
merical cells is set equal to the local fluid velocity. 
TESS is capable of solving general systems of hyper- 
bolic equations in conservation-law form, in one, two, 
and three dimensions, with arbitrary coordinate geom- 
etry. TESS is also a parallel code, configured to run 
on distributed-memory supercomputers, achieving paral- 
lelism via domain decomposition. It is therefore capable 
of very high resolution calculations. TESS is ideal for 
the study of fiuid disks whose bulk motion is supersonic, 
especially for problems which need high accuracy which 
could potentially be compromised by large advection er- 
rors. For more detai l s on t he numerical method, see 
iDuffell fc MacFadven I ([201 1[ ). All of the features most 
recently added to the code (e.g. 3D and parallelization) , 
will be detailed in a future publication. 

Ordinarily, TESS accomplishes its mesh motion by cre- 
ating a Voronoi tessellation of the computational do- 
main to determine the shape and size of its finite vol- 
umes. However, the TESS scheme is modular so that 
tessellations particularly suited to a given flow can be 



4 



Duffell & MacFadyen 



implemented. For calculations involving proto-planetary 
disks, we choose to divide the computational domain into 
annular segments (similar to the tessellation which is per- 
formed by standard codes using polar coordinates) and 
allow the annular segments to each independently move 
with the local fluid angular velocity (Figure [2|) . Since 
the vast majority of the fluid motion is Keplerian, this 
scheme is very close to a Lagrangian formulation; there 
will be some advective fluxes in the radial direction, nec- 
essarily, but in principle we can always arrange, for ex- 
ample, to have the faces move in such a way that there 
is exactly zero advective flux in the azimuthal direction. 
In the present work, we simply give each cell the local 
Keplerian orbital velocity. This choice removes the bulk 
background flow and allows for computation of the fluid 
flow to take place on a grid effectively co-moving with 
the fluid. As such, subtle flow features are captured 
which would be artiflcially diffused and dissipated if a 
flxed mesh were used. 

The numerics are essentiall y identical to the TESS code 
puffeh k MacFadven1l20ll[ ): only the shape of the cells 
has changed. The reason for this specialization in the 
case of disks is that while TESS is excellent at reducing 
diffusion due to large bulk motions, it can suffer from a 
small amount of numerical noise in the presence of large 
shearing motion; when two cells shear past one another, 
the face shared between the cells rotates quickly, produc- 
ing noise. While this noise is a very small price to pay for 
the vast reduction of diffusive fluxes, it is best to avoid it 
if possible. The disk-tessellation into annular segments 
we have chosen for this study is an excellent solution to 
this problem. Also, because it is a more restrictive ge- 
ometry, the "tessellation algorithm" is extremely efhcient 
and straightforward to implement, when compared to a 
Voronoi tessellation. We note that this idea is similar in 
spirit to orbital advection schemes like FARGO (Masset 
[200Qi) . but that the numerical method is clearly distinct. 
FARGO is designed to work with the more restrictive log- 
ically cartesian grids, whereas our zones have a variable 
number of neighbors, usually about six in 2D. 

All the calculations were done using the basic TESS al- 
gorithm, but there are a few differences because we were 
able to take advantage of the simpler geometry. The 
piecewise linear reconstruction has been modifled for the 
annular segments, so that it looks exactly like a standard 
piecewise-linear method in the (j) direction. The other 
important distinction is that we formulate the hydrody- 
namic equations in terms of the angular momentum, so 
that angular momentum is explicitly conserved, as is de- 
sirable for disk calculations. 

3.2. Initial Conditions 

All of our calculations are specifled completely by the 
planetary mass as a fraction of the gap threshold mass 
Mp/MTh and the disk Mach number A4, which in this 
work will always be set equal to 20. This value is con- 
sistent with observations of proto-planetary disks. The 
initial conditions are given by: 



CpSo(r, 0)/7 



(21) 



Cp — Vtprp/M. 



(18) 

(19) 
(20) 



For simplicity, we choose a uniform density and pres- 
sure proflle, though arbitrary profiles can be imple- 
mented straightforwardly. Also for simplicity, we leave 
the planet at a fixed radius, though migration can eas- 
ily be taken into account. The exact value of Ep is 
arbitrary, because the disk is not a source of gravity 
in our calculations; a change in Ep is equivalent to a 
change of units. Note also that our initial conditions 
have included a planetary atmosphere (an overdensity 
corresponding to the exponential in the formula for the 
surface density; $p(r, 0) is the planetary potential, de- 
scribed in the next section). This choice was not only 
convenient for reducing transients, but actually neces- 
sar y for capturing the dyna r nics co rrectly. It was noted 
bv jPong. Rafikov fc Stone I (120111: ) that orbital advec- 



tion schemes like FARGO (iMasset 



20001) can potentially 



arrive at misleading results because these codes attempt 
to "cheat" the Courant condition for their time-steps 



(22) 



by setting Xmax to the sound speed rather than the or- 
bital velocity (this is one of the primary motivations of 
FARGO). The problem cited by Dong et al. is that there 
is also a gravitational timescale associated with accre- 
tion onto the perturbing planet, so that if a code's time- 
step criterion fails to take into account this additional 
timescale, qualitatively incorrect dynamics can result 
(e.g. spurious gap opening). We have found this to not 
pose a problem for our numerical method, if we include 
the planetary atmosphere as an initial condition. This 
is because we are searching for the quasi-steady-state so- 
lution, which should not have any features which evolve 
on the gravitational timescale. Because we choose initial 
conditions including the planet's atmosphere which are 
close to stationary, the system can be evolved using the 
much longer timesteps associated with the sound crossing 
time of a computational zone. Note that this improves 
computational efficiency by at least a factor of the Mach 
number, which in this work is 20. In fact, because the 
zones at the inner boundary have a much larger orbital 
velocity, the speed-up is actually about twice this factor. 

3.3. Planetary Potential 

Because the planet's position is inside the computa- 
tional domain, it is not possible to use the exact (diver- 
gent) potential 



(23) 



where s is the distance from the planet. Instead, it is 
standard practice to use an approximate potential with 
a smoothing length, eg. The potential should quickly 
converge to 1/r at distances larger than Cg, and so we 
use a potential which converges at nth order: 



$(")(s) 



(s" + ej)V" ' 



(24) 



where s — \ fp\ is the distance f rom the fluid element 
to the planet. Typically we follow iDoiiget al. I (|2011aD 
and use the fourth order version of this potential, n = 4 
and choose Cs/rp = .005 (1/10 of a scale height), which 
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Fig. 3. — Position of the spiral density wave in the disk. The 
curve plotted is the analytic prediction l[9jl which assumes a linear 
perturbation. Points are density maxima in the numerical calcu- 
lation. The analytic formula is only valid far from the planet. In 
the right panel, we show a close-up of the planet. Here we see the 
same deviation from this curve as did iDong et al71 Il2011al) (Paper 
I, Figure 2). 

is sufficient to produce the correct density waveform, as 
shown in section 23] 

3.4. Boundary Conditions 

Boundary conditions are particularly challenging for 
astrophysical disk calculations. There necessarily exists 
an inner boundary, at which the gravitational force due 
to the parent star is strongest, and numerical integrations 
are most prone to inaccurate cancellation between cen- 
trifugal and gravitational forces. In general, if the system 
is allowed to evolve for long enough, spurious errors will 
develop at this inner boundary. If the boundary condi- 
tion there is not well behaved, these errors will eventually 
propagate to fill the rest of the computational domain. 
Ideally, we would like a boundary condition that allows 
waves to simply propagate out of the inner boundary 
without refiecting back. Unfortunately such a boundary 
condition is ill-posed. 

Rather than trying to directly implement a compli- 
cated inner boundary condition, we found an effective 
way to prevent reflections by exponentially damping per- 
turbations within some given radius r < Tmin- Thus, 
waves propagating toward the inner boundary are com- 
pletely damped and therefore unable to reflect back. 
f-min behaves as the boundary of our computational 
domain, and we interpret this damping as an outflow 
boundary condition at this radius. As we shall see, 
the relevant dynamics are not negatively impacted by 
the simpliflcations we make at radii r < rmin- Due to 
the good behavior of this boundary treatment, we de- 
cided to use the same boundary condition for large radii, 
r > Tmax- Typically we use the values r,„i„ = QAr^, 
fmax = 3rp, whereas our total computational domain ex- 
tends to 0.25rp < r < Avp. 

4. RESULTS 

iDong et al. I H2011aD warn that this is a challenging 
problem which requires very high resolution and solid 
numerics in order to get it right. The reason for this 
is that everything hinges upon calculating the proflle 
of the linear density wave very accurately, and cap- 
turing the extremely weak shock which subsequently 
forms. In fact, even to get the calculation right in 
the linear approximation using a discrete fourier trans- 
form, .GoMmaniLRafikmL, (,2001i) required 4096 x 8192 
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Fig. 4. — Density wave at t = 1.89, for planet mass 
-A/p = .0209-Mt(,, co mpa r ed wi th the semi-analytical predictions 
of Go odman fc RafikoTI ||200D. The sohd curve is Goodman & 
Rafikov's calculation, and the dashed curve corresponds to our nu- 
merical meas urement. Variables r) and x given by equations 
mul l and mil l, respectively, x is proportional to the difference in 
surface density from the background value, and r) is related to the 
angular position. The amplitude of our wave was reduced by 8% 
in this plot to account for a discrepancy in normalization. 

points in k-space. Using the Athena code, iDong et al. I 
()2011al) were able to capture the dynamics accurately 
in the shearing box approximation, which required 3072 
X 16384 grid cells. We are performing a global calcu- 
lation, which will require even higher resolution, if only 
because our computational domain is significantly larger. 
This fact alone should increase the computational cost by 
roughly a factor of five in resolution. Typically we use 
14400 cells in the radial direction, and roughly 27r times 
this number in the azimuthal dimension. Note that the 
number of azimuthal cells varies with radius, which is an- 
other feature of our method. We choose our azimuthal 
resolution so as to keep a fixed cell aspect ratio of 1:1. 
The total computational domain consists of roughly 500 
million zones. 

Because this is a challenging problem, we first demon- 
strate that our method accurately captures the predicted 
shock formati on, essentially reproducing the s hearing- 
box results of IDong. Rafikov fc Stone] (|2011b[) , but in 
the global case. 

4.1. Linear Phase and Shock Formation 

Equation ([9]) predicts the position of the spiral den- 
sity wave, which we compare with our numerical cal- 
culation in Figure |3l To compare the analytic formula 
with our results, we measure the value of (j) at which the 
surface density peaks for a given radius. This is a rela- 
tively easy test, for which we found good agreement even 
at low resolution (for example. Figure [Tj was computed 
at lower resolution, using only 4096 radial zones). Very 
close to the planet, our result deviates fro m this predic- 
tion, in the same respect that was found bv IDong et al. I 
(|2011al) . Actually, the density maximum is not predicted 
to coincide with this curve exactly; this would correspond 
to the density wave peaking at 77 = 0, which is not pre- 
dicted to occur by any theory. However, we don't have to 
move far from the planet before the value of 77 at which 
the density peaks is much smaller than the angular posi- 
tion |77mox| << M.cj), meaning the position of the density 
maximum is generally a good proxy for "wave position" . 
Regardless, the deviation close to the planet is expected 
and was also seen in the shearing box calculations. 

Before nonlinear evolution sets in, the shape of 
the wave should agree with the linear calculations by 
iGoodman fc Rafikov I (|200lD . because the shearing box 
approximation is still appropriate in this regime. Our 
calculation of the density wave agrees with linear the- 
ory (Figure |4]) up to an overall normalization, which is 
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Fig. 5. — The nonlinear evolution of the density wave into a shock (Mp = .0209MTh)- The left and center panels show th e tw o sep arat e 
shocks that form (inner disk and outer disk, respectively). We plot the perturbation using the variables rj and Xi equations l|10| l and jlTJl, 
resp ectively. The shock is shown at radii corresponding to i = 1.89, 10.8, 29.7, 60.8, 106, and 246. The function i(r) is given by equation 
111211. The right panel shows the scaling relations for the width and the height of the shock. This figure can be compared with Figure 1 of 
IDong. Hafikov Ik Stone"! Il2011bl ) (Paper II), where they performed the same calculation in the shearing box approximation. 
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Fig. 6. — Normalized torque density for a planet of mass 
Mp = .0209Mj'h- Normalized torque density is defined as 
dT/dr X Mxh/i-M'^oGMp).. The solid curve is a linearized cal- 
culation by Rafikov in a shearing-box. Dashed curves are our nu- 
merical result for the inner and outer disk. 



arbitrary and difficult to determine exactly when tak- 
ing into account the planetary atmosphere introduced in 
our initial conditions (j2ip . The shape of the waveform 
is correct, though the normalization of the amplitude is 
8% greater than shearing box predictions. 

The most important test of our calculation is the ac- 
curate capturing of the shock, since this is the mech- 
anism for exchange of angular momentum between the 
planet and the disk. In figure [5] we track the shock for- 
mation. Our results ar e consis tent with the results of 
IDong. Rafikov fc Stone 1 (|2011b( ) for this calculation. We 
also confirm the scaling relations predicted by Burgers' 
equation for the amplitude and width of the dissipating 
shock as a function of t(r) (Figure [SI right panel). 

4.2. Torque Density 

Before reporting the results of the fully nonlinear cal- 
culation, we discuss the torque density in the linear per- 
turbation. Planetary migration timescales can be calcu- 
lated by adding up the total torque on the planet due to 
the net gravitational pull of the perturbed disk density 
profile. In analytic calculations, this is typically done in 
Fourier space, but in this numerical study we can directly 



calculate torque density as a function of radius: 

^ = Js^rUrdcf^ (25) 

where is the (j) component of force per unit mass be- 
tween the planet and the disk. 

Early analytic results such as iGoldreich &: Tremaine I 
(flQSQ) calculated the asymptotic behavior of the torque 
density, and found it to decay like a power-law. Re- 
cently, however, this asymptotic behavior was shown by 
IDong et al. I ()2011aD to be incorrect. They found that 
the torque density changes sign at a particular distance 
from the planet, 

r_ = rp ± 3.2/1 (26) 
which corresponds to 

= 17. (27) 

iRafikov fc Petrovich I (|2012| ) explained this disagreement 
by doing a more careful analytical calculation, without 
assuming that different Fourier harmonics were confined 
to the vicinity of Lindblad resonances. Figure [5] shows 
the torque density in our global calculations. As in other 
recent works, we observe a change in sign of torque den- 
sity at finite distances from the planet, |r— r^j — 2.8ft. and 
3.6/i, for the inner and outer disks, respectively. Both of 
these radii are in agreement with the analytic value of 
t- = 17 found in the shearing box case. Actually, this 
particular result was not derived for the global case, but 
we find that expressing r_ in terms of the t coordinate 
gives a reasonable prediction for where the torque den- 
sity becomes negative for both the inner and outer disk 
in the global case. 

The difference between the torque in the inner and 
outer disk gives the net Lindblad torque on the planet. 
Since we have a converged global calculation of torque 
density, we can numerically integrate it to find the net 
torque. We find the net torque to be 

T = -0.1 r^^aGMl/MTh. (28) 

Note that in reporting this final result we have removed 
the scaling with Mach number, because this calculation 
was specific to the case M. = 20. Additionally, this result 
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Fig. 7. — Normalized angul ar m omentum flux ^{r)/^(rp) for Mp 
the predicted scaling relation jlTJ. 



.0209Mt,i and Mp = .24AfTh(lOM0). Thick dashed curves show 




Fig. 8. — A IOMq planet giving hints of gap formation. Color 
represents the perturbation to the surface density, (5E/So. Dotted 
line s ar e the theoretically predicted radii for shock formation, given 
by I I13II . A downsampled subsection of the computational domain 
is showrn, .75rp < r < 1.25rp. 

assumes uniform density and pressure in the disk, and 
would be modified significantly in the presence of non- 
trivial density and pressure profiles. Because the planet 
is at a fixed radius, we are able to conpute only the Lind- 
blad torque, and not the corotation torque. 

4.3. Angular Momentum Flux 

From our global numerical calculations we can con- 
firm the prediction of g lobal angular m omentum flux 



lirm tne prediction ol g lobal angular m omentum nux 
as a function of radius ( Rafikov 11200231 ) . In Figure [7] 
we show the function ([H]) as computed from our data. 
For comparison, w e show the theor etical scaling relation 
pT|) calculated bvl Rafikovl ()2002aD . We appear to have 
better agreement with the semi-analytic theory in the 



outer disk than the inner disk, but the overall picture 
is clear; the angular momentum flux is roughly uniform 
between the planet radius and the radius at which the 
shock forms, after which shock dissipation causes the per- 
turbation to deposit its angular momentum. It should be 
noted that since pTl) is a scaling relation, we could get 
the initial waveform completely wrong and still have the 
correct (normalized) angular momentum flux, as long as 
the shock is correctly captured. Predictions involving the 
angular momentum flux are therefore robust. 

4.4. Gap Formation 

Because the flux of angular momentum is not uniform, 
we see time-dependent disk evolution, which can lead to 
gap formation. In Figure |8l we show the perturbation 
due to an intermediate- mass planet, .24Mt;i = lOAf^, 
after 100 orbits. Two gaps appear to be forming, one 
at the inner shock position and one at the outer shock 
position. The question arises as to how long the disk 
evolution should continue; this is a question we cannot 
answer directly, but because the gap opening is due to 
shock dissipation, there is no reason to expect this trend 
to saturate anytime soon. Because our "gap" is only a 5% 
deviation in the surface density, it is not clear what this 
system will look like at late times, but after 100 orbits 
the surface density in this region still has a negative time 
derivative, showing no signs of a plateau. The time-scale 
for gap opening is measured to be roughly 1500 orbits. 
Regardless of what happens at late times, we have cer- 
tainly found that a lOM^ planet is massive enough to 
drive material away from it on a 1000 year timescale. 
However, since we have only watched it for 100 years, 
it is not immediately obvious that this trend continues 
for 1000 years, though with negligible viscosity it is the- 
oretically predicted that th is tren d will persist until the 
gap is very deep ( RafikoTI [2002bl ) . IMuto et al. 1 H2010f ) 
saw a similar pattern in the shearing-box; their calcu- 
lations suggested that the density perturbation became 
nearly statio nary after a f ew hun dred orbits, but it was 
suggested bv lDong et al. I ()2011af ) that these calculations 
could have been underresolved. 

We should note that this gap formation is the result of 
an accumulation of angular momentum flux over many 
orbits. For this reason , it wo uld not be possible for 
iDong. Rafikov fc Stone I ()2011bD to have observed gap 
formation in their calculations, because they used an in- 
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flow boundary condition in their shearing box. If they 
had used a periodic boundary condition, and if they had 
attempted using a lOM^ planet, we expect that they 
would have seen the same gap that we have found in the 
global case. 

5. SUMMARY 

We have performed high accuracy global protoplan- 
etary disk c alculat ions which extend the results of 
iDong et al. I (|2011al ) and have confirmed the analyti- 
cal predictions of Rafikov for the case of a global disk. 
Specifically, we directly measured the linear waveform, 
shock formation, and ang u lar mo mentum flux, and find 
agreement with iRafikov 1 ()2002al ). We have also cal- 
culated the torque density for a low mass perturber 
and found agreemen t with recent shearing-box r esults 
(|Dong et al. I l2011al: IRafikov fc Petrovich I [20121 ) that 
find torque density to become negative a few scale heights 
away from the planet. 

We have shown directly that a perturbing mass 
Mp < Mrh (specifically, ten earth masses) can have a 
non-trivial impact on its environment, potentially lead- 
ing to gap formation. All of this, of course, assumes that 
the thin-disk approximation is valid to use here, which 
may not be true close to the planet, but it may be a rea- 
sonable assumption where the shock forms. How much 
this assumption affects the result is unknown; this should 
be checked by three dimensional simulations. We have 
avoided speaking specifically about what conditions are 



necessary to open a gap. To directly explore gap-opening 
criteria, we need to include planet migration and disk vis- 
cosity in our calculations, which is the topic for future 
work. 

We have also demonstrated the power of using a mov- 
ing mesh for performing hydrodynamical calculations in- 
volving gaseous disks, especially when a cylindrical grid 
geometry is used. This result is the first of many plane- 
tary migration calculations using this method. Now that 
we have demonstrated that it is possible to capture all 
of the important and subtle details of the shock forma- 
tion and we have presented a converged global result, we 
can begin to study planet migration in this intermediate- 
mass range. In addition, it is clear that this method 
should prove itself useful for studying many other sub- 
jects, for example binary black hole systems, the magne- 
torotational instability, and accretion onto compact ob- 
jects. 
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